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Abstract 

Vapor-liquid menisci of the truncated and shifted Lennard-Jones fluid between parallel 
planar walls are investigated by molecular dynamics simulation. Thereby, the characteristic 
energy of the unlike dispersive interaction between fluid molecules and wall atoms is system- 
atically varied to determine its influence on the contact angle. The temperature is varied as 
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well, covering most of the range between the triple point temperature and the critical temper- 
ature of the bulk fluid. The transition between obtuse and acute angles is found to occur at a 
temperature-independent magnitude of the fluid-wall dispersive interaction energy. On the ba- 
sis of the present simulation results, fluid-wall interaction potentials can be adjusted to contact 
angle measurements. 

Introduction 

A major challenge for molecular modeling consists in the definition of unlike interaction potentials. 
In the past, a variety of combination rules were proposed, none of which was found to be valid in 
general. Several of these, including the Lorentz-Berthelot combination rule, are considered to 
be a good starting point for further adjustment in most cases'^. The present work contributes to 
understanding the dispersive interaction between a solid wall and a fluid, which is essential for the 
analysis of adsorption and microscopic flow properties. 

In principle, the Lorentz-Berthelot rule can be applied for effective pair potentials acting be- 
tween the fluid particles and the atoms of a solid walll^El, based on size and energy parameters 
derived from properties of the solid and the fluid. However, while using combination rules to ex- 
trapolate from homogeneous bulk solid and fluid properties to interfacial phenomena can lead to a 
good agreement with the actual behavioi^, this approach has only shaky theoretical foundations'^. 
Unlike pair potentials between a fluid and a solid can only be expected to give reliable results 
if they are developed using actual information on fluid-wall contact effects. Since adsorption in 
nanopores can be studied on the basis of effective pair potentials'^'^, it is obversely possible to 
fit model parameters to adsorption isotherms'^. The present study follows the line of research, 
suggested by Werder et al^, of adjusting unlike parameters to contact angle measurements. 

Using density functional theory, Teletzke et al^ examined the dependence of wetting and 
drying transitions on characteristic size and energy parameters of the fluid-wall dispersive inter- 
action. Subsequently, SOKOLOWSKI and Fischer'^ as weU as Giovambattista et a/.^^^ in- 
vestigated fluid density profiles in extremely narrow channels for several values of the fluid-wall 
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dispersive energy and surface polarity, respectively. On the microscopic and the nanoscopic level, 
the statics and dynamics of fluids under confinement and the corresponding three-phase contact 
lines can also be investigated by the lattice Boltzmann method™. 

Molecular dynamics (MD) simulation can be applied to this problem as well, leading to a con- 
sistent molecular approach. The increase in computing power and the development of massively 
parallel MD software provide tools for simulating systems with a large number of particles. Sys- 
tem sizes accessible to MD simulation are getting closer to the smallest experimental settings. This 
allows comparing simulation data directly to experimental observables for a growing spectrum of 
properties, including the contact angle. The truncated and shifted Lennard-Jones (LJTS) poten- 
tiaP is used in the present work for describing both the fluid-fluid and the fluid-wall interaction, 
leading to systems that extend previous studies on interface properties for the LJTS fluid™. 
Hysteresis ^-^ as well as the formation and growth of liquid precursor layers on the surface^^are not 
discussed as the present work deals with equilibrium properties of the phase boundary only. 



Simulation method 

Like the original Lennard-Jones (LJ) potential u^\r) = 4e [(a/r)^^ - (a/r)^] , the LJTS model^^ 



, I u^\nj)-u^\r,) rij<r, 

Uij{rij) ={ (1) 



Tij > rc, 



with a cutoff radius of rc = 2.5 cr, accurately reproduces the thermophysical properties of several 
non-polar fluids, in particular noble gases and methane, when adequate values for the size and 
energy parameters o and e are specified^*^. Due to the relatively small cutoff radius, molecular 
simulation is comparably fast, while the full descriptive power of the LJ potential is retained even 
for systems with phase boundaries'^. 

In order to accurately describe properties of a solid material, it is usually necessary to employ 
multi-body potentials which have a large number of model parameters and are computationally 
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quite expensive " I The present study, however, does not regard the properties of a specific wall 
material but rather the influence of the fluid-wall dispersive interaction on the fluid itself. Ac- 
cordingly, a layered wall was represented here by a comparably straightforward system of coupled 
harmonic oscillators, using different spring constants Dy and D^^ for the transverse vibration with 
respect to the layers 

u]{yd = ^{y-y-)\ (2) 

wherein y° is the equilibrium value of the y coordinate (i.e. the direction perpendicular to the wall), 
and longitudinal oscillations, 

u^j{rij) = ^{nj-A)\ (3) 

with respect to the equilibrium bond length A between neighboring atoms i and j. 

Fluid-wall interactions can be represented by fulP or slightly modified LJ potentials, acting 
between fluid particles and the atoms of the solid. Following this approach, the LJTS potential 
with the size and energy parameters afw = C7 as well as 

Cfw = Ce, (4) 

was applied for the unlike interaction using the same cutoff radius as for the fluid. Potential pa- 
rameters for the molecular models of the fluid as well as the solid component were chosen such as 
to represent methane and graphite, respectively. For the fluid, the LJTS size and energy parameters 
o = 3.7241 A and £/k= 175.06 K, as well as the molecular mass m = 16.04 u were usedl^, so that 
the carbon-carbon (C-C) bond length in graphite A= 1.421 A^^can be expressed in LJ units as A 
= 0.3816 o, while the interlayer distance Y = 3.35 A corresponds to 0.8996 o. 

The spring constant Dxz = 15600 N/m related to the sp^ bonds was adjusted to the C-C radial 
distribution function obtained from simulations with the Tersoff'^ potential. This distribution 
had to be rescaled because as previously shown'^^'^^', the Tersoff potential deviates by about 3% 
from the actual bond length in graphite. In agreement with the relation between the C-C bond 
energy (4.3 eV) and the interaction energy between adjacent graphite layers (0.07 eV)'^, the in- 



Figure 1: Simulation snapshots for the reduced fluid-wall dispersive energy C, = 0.09 (left) and 
0.16 (right) at a temperature T = 0.73 e/k. The upper half is reproduced in the bottom to illustrate 
the effect of the periodic boundary condition. 

terlayer spring constant was specified as Dy = Dxz/ 60. Equations ([2j) and ^ ensure that the wall 
atoms oscillate around the fixed y coordinate that corresponds to their layer, while no particu- 
lar X and z coordinates are preferred because the atoms are only connected with their immediate 
neighbors to permit individual sliding of the wall layers. 

Massively parallel canonical ensemble MD simulations were conducted with the program isl 
Mardyr^ to obtain the contact angle dependence on the temperature and the reduced fluid-wall 
dispersive energy. For all simulation runs, the integration time step was specified as 1 fs. Vapor 
and liquid were independently equilibrated in homogeneous simulations for 10 ps. This was fol- 
lowed by 200 ps of equilibration for the combined system, i.e. a liquid meniscus surrounded by 



vapor, with a wall consisting of four to seven layers, cf. [Figure 1[ where the starting configura- 
tion contained a planar vapor- liquid interface perpendicular to the z coordinate axis. Note that the 
distance from the wall is given by the y coordinate, while z is the characteristic direction for the 
density gradient of the fluid. The periodic boundary condition was applied to the system, leaving 
a channel for the fluid with a height of 27 o between the wall and its periodic image. 

Via binning, the density profiles were averaged over at least 800 ps after equilibration. The 
arithmetic mean density pvi = (p' + p")/2 was applied to define the position of the phase boundary, 
where p' and p" are the saturated bulk densities of liquid and vapor which are known for the LJTS 
fluid from previous work^^. In the immediate vicinity of the wall, the fluid is affected by short- 
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range ordering The influence of this phenomenon was minimized by taking density 

averages over a bin size of about 1 a, cf. Figure 2 following GlOVAMBATTlSTA et al. . A circle 



was adjusted to the positions of the interface in the bins corresponding to distances between 3 and 
11a from the wall, and the tangent to this circle at a distance of 1 a from the wall was consistently 
used to determine the contact angle. 
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Figure 2: Vapor-liquid interface profiles for the reduced fluid-wall dispersive energy C, = 0.07 (up- 
ward triangles), 0.10 (squares), 0.13 (circles), and 0.16 (downward triangles) at a temperature T 
= 0.82 e/k. Note that the full lines are circle segments, adjusted to the data points that are repre- 
sented by symbols. The almost perfect match between the individual simulation results, indicated 
by their collective agreement with the fit, reflects the precision of the present simulation data. At 
T = 1 e/k, however, the margin of error becomes more significant. 



Simulation results 

Menisci between parallel planar walls were simulated for a reduced fluid-wall dispersive energy C, 
between 0.07 and 0.16 at temperatures of 0.73, 0.82, 0.88, and 1 e/k. Note that for the bulk UTS 
fluid, the triple point temperature is about 0.65 e/k according to VAN Meel et al.^^ while the 
critical temperature is 1.0779 e/k according to Vrabec et al.^^, so that almost the entire regime 
of vapor-liquid coexistence was covered here. 

High values of ^ correspond to a strong attraction between fluid particles and wall atoms, 
leading to a contact angle < 90°, i.e. to partial (t? > 0°) or perfect (t} = 0°) wetting of the 

6 



surface. As expected, with increasing fluid-wall dispersive energy, the extent of wetting grows, 
cf. Figure 1 As can be seen in Table 1 and Figure 2\ the transition from obtuse to acute contact 
angles occurs at C, values between 0. 1 1 and 0. 13 over the whole studied temperature range. Present 
simulation results were correlated by 



cosi?(r,C) 



1 + 



.1.7 



27 



tanh 



0.087' 



(5) 



where t = {I —T /Tc)~^ approaches infinity for T — > Tc, while Z = 0.119 is the reduced fluid-wall 
dispersive energy that leads to a contact angle of = 90°. 

Perfect wetting or drying is present where Eq. ^ yields cos > 1 or cos < — 1 , respectively. 
In particular, both the value of Z and the symmetry relation 



cos-i?(r,Z-AO = -cos-i>(r,Z + AC), (6) 
were found to be temperature-independent. Note that the simulation results stronlgy suggest such 



a symmetry property, cf. Figure 3[ which is not an artifact of the correlation. 



Figure 3 shows, particularly for high temperatures, that there is a narrow range of ^ values 



that lead to the formation of a contact angle as opposed to perfect wetting or drying. The present 

Table 1 : Contact angle of the LJTS fluid on graphite from MD simulation as a function of reduced 
fluid-wall dispersive energy and temperature. 
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Figure 3: Simulation results and correlation, cf. Eq. (|5J), for the contact angle in dependence of the 
reduced fluid-wall dispersive energy C, at temperatures T = 0.73 (diamonds and solid line), 0.88 
(squares and dashed line) as well as 1 e/k (circles and dotted line). 
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Figure 4: Simulation results and correlation, cf. Eq. ([5|, for the contact angle in dependence of 
the temperature at fluid- wall dispersive energies ^ = 0.09 (upward triangles and solid line), 0.10 
(squares and dashed line), 0.12 (bullets and solid hne) as well as 0.14 (downward triangles and 
dotted line). The entire temperature range between triple point and critical point of the bulk fluid 
is shown. 



plots agree qualitatively with those determined by Giovambattista et a/.^^for the influence of 
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the polarity of hydroxylated silica surfaces on the contact angle formed with water. In Figure 4 



can be seen that the extent of wetting (for ^ > Z) or drying (for C, < Z), respectively, increases as 
the temperature approaches Tc. Eventually, this leads to the known phenomenon of critical point 
wetting!^ for the whole range above a wetting temperature T^. The effects described above can be 
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accounted for by the Young equation 

cos.(r.C) = '^''"-^-P";-f"-^-P'> . (7) 

which relates the vapor-liquid surface tension 7 and the contact angle to the interfacial tension 7s 
that acts between wall and vapor or liquid. The deviation from = 90° increases with T 
because y{T) converges to zero faster than the density difference p'{T) — p"{T) between the two 
fluid phases. The relevant critical exponents in case of the LJTS fluid are 

,._^in[p'(r)-p"(r)]^l 

T^T, in(rc-r) 3 

for the saturated densities of the bulk fluidl^, in accordance with the Guggenheim'^ approach, 
and 

-1.21, (9) 



d\n{T^-T) 

for the vapor- liquid surface tension!^, confirming a similar value (1.26) obtained from fluctuation 
theory of critical phenomena'^. 

The correlation given by Eq. ([5]) suggests for the present system a first-order transition between 
partial and perfect wetting or drying, respectively, as described by CahnES. With C -Z, |Fig 



ure 3 shows that the contact angle ■& depends linearly on the fluid-wall dispersive energy, and the 



symmetry property suggests for ^ = Z that 7s does not depend on the density of the fluid. 



Comparative discussion 

Qualitatively, the symmetry relation given by Eq. ^ corroborates MonsonI^ who obtained the 
same property based on mean-field DFT calculations. If each layer of the wall is approximated as 
a plane of uniform density r] , the well depth of the fluid- wall dispersive interaction, which can be 
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used to compare different interaction models quantitatively, is given by 



(10) 



for a system of L layers with an interlay er distance of F. In the present case with rc = 2.5 cr for the 
LJTS potential as well as a surface density of r/ = 5.287 a^^, an interlayer distance of 7 = 0.8996 
a, and the number of layers L = 3 (the fourth layer of the wall is beyond the cutoff radius), one 
obtains 

= 17.29 ytTc ■ (11) 

normalized by the critical temperature Tc = 1.0779 e/fc of the bulk LJTS fluidl^. The transition 
from obtuse to acute contact angles, occurring at Z = 0.119 in the present case, therefore corre- 
sponds to a well depth of Wz = 2.057 fcTc. 

A very similar system was investigated by BUCIOR et alW, albeit for the LJTS fluid with a 
cutoff radius of 2^/^ a as well as a rigid wall with a surface density ^ = \ .7342 <5^ and a single 
solid layer, i.e. L = 1. For that fluid, the critical temperature is given by SCHRADER et al^^ as 
Tc = 0.9999 ejk, so that the fluid-wall dispersive energy was related to the well depth by = 



5.429 fcTc ■ C, according to Eq. (10). From the density profiles of BuciOR et al. one finds that 
a rectangular contact angle is reached for a reduced fluid- wall dispersive energy C, between 0.61 
and 0.7, i.e. 3.3 kT^ <Wz < 3.SkTc which is on the same order of magnitude as the present result. 
The quantitative deviation has to be attributed to the different solid structure, since the single wall 
layer of BUCIOR et al.^ leads to a faster decay of the fluid-wall dispersion with respect to the 
distance than in the present simulations where three layers can directly interact with the fluid. A 
larger value of Wz/{kTc) is required to compensate for the effectively smaller length scale of the 
dispersive interaction. 
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The present correlation, cf. Eq. ([5]), predicts perfect wetting for 



r„(i|l= I 1-0,144 



tanh 



W 

QMS- 1.37 



1 

- 1 




(12) 



The transition to perfect wetting was also simulated by Bo JAN et al.^^ who applied the Monte 
Carlo method in the grand canonical ensemble to neon on metal surfaces. Thereby, the full LJ 
fluid, with Tc = 1.310 e/fc as determined by LOTFI et al.^^, was used to model neon. With W 
= 2.13 kTc, representing magnesium, they obtained a wetting temperature of Tw ~ 0.50 Tc, as 
opposed to the present results which imply that perfect wetting is only reached above 0.974 T^, in 
the immediate vicinity of the critical temperature. 

As BOJAN et a/. 1^ themselves remark, their calculations predict a much lower wetting tem- 
perature than a similar previous study by SOKOLOWSKI and Fischer'^ on the local structure of 
fluid argon in contact with solid carbon dioxide (with the same value of W /e). The latter MD 
simulation results used a size parameter C7fw/<7 that is significantly smaller than unity in case of 
argon and carbon dioxide'^^'^and therefore cannot be directly compared to the present work, since 
varying the interaction length scale can lead to qualitatively different properties such as a change 
in the order of the wetting transition!^. However, it should be pointed out that significantly better 
agreement was obtained here with the results of SOKOLOWSKi and FlSCHERl^than with those of 
Bo JAN et al.^. 

The results presented above can now be used to provide an estimate for the magnitude of the 
dispersive interaction between fluids and wall materials for which experimental data on the contact 
angle are available. For instance, regarding the refrigerant R134a (1,1,1,2-tetrafluoroethane) at 
temperatures between 10 and 80 °C, Vadgama and Harris"^^ obtained contact angles of 5.5° ± 
r on copper and T ± 1° on aluminum. With T = 0.85 T^, which is 45 °C for R134a, the weU 
depth can be estimated as W /k ^ 2.9 Tc= 1 100 K in both cases on the basis of Eqs. (|5| and (11). 
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This can be related to a molecular model of the dispersive interaction by 



W^ = -p^min/ dXlnX f dv Ui^(\{y + vf + X^V'^\ , (13) 
>'>OJo Jo ^ / 

wherein Pw is the density of the solid wall and Mf^ is the dispersive interaction potential acting 
between a fluid molecule and a wall atom. Note that, while Eq. ([10]) corresponds to a sum over 



truncated and shifted LJ-10-4 terms, Eq. (13) does not rely on any particular assumption on the 
internal structure of the solid wall. If e.g. a LJ-12-6 potential is used for Wf^, it corresponds to 
a LJ-9-3 interaction. This reasoning can plausibly be applied to all fluids that do not exhibit an 
excessively polar or anisotropic structure. 



Conclusion 

The contact angle formed between a wall and a vapor-liquid interface was determined by canonical 
ensemble MD simulation, where the magnitude of the dispersive fluid-wall interaction and the 
temperature were varied. Over the whole temperature range under investigation, the contact angle 
dependence on the fluid-wall dispersive energy was found to follow a simple symmetry law. At a 
temperature-independent value of the reduced fluid-wall dispersive energy, the interfacial tension 
between vapor and solid as well as hquid and sohd is equal, corresponding to the transition between 
obtuse and acute contact angles. 
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